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Our recently developed variant of variationnally optimized perturbation (OPT), in particular 
consistently incorporating renormalization group properties (RGOPT), is adapted to the calculation 
of the QCD spectral density of the Dirac operator and the related chiral quark condensate (qq) in 
the chiral limit, for nf = 2 and n/ = 3 massless quarks. The results of successive sequences 
of approximations at two-, three-, and four-loop orders of this modified perturbation, exhibit a 
remarkable stability. We obtain (qq)l/ 3 =2 (2 GeV) = —(0.833 — 0.845)A2, and (qq)l/^ =3 (2 GeV) = 

— (0.814 — 0.838)A 3 where the range spanned by the first and second numbers (respectively four- 
and three-loop order results) defines our theoretical error, and A n/ is the basic QCD scale in the 
MS- scheme. We obtain a moderate suppression of the chiral condensate when going from nf = 2 to 
n/ = 3. We compare these results with some other recent determinations from other nonperturbative 
methods (mainly lattice and spectral sum rules). 
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I. INTRODUCTION 

The chiral quark condensate ( qq) plays a central role in QCD nonperturbative dynamics, being together with 
the pion decay constant the other principal (lowest dimensional) order parameter of spontaneous chiral symmetry 
breaking, SU(rif)i, x SU(nf)n —> SU{nf)v for rif massless quarks, the physically relevant cases being n/ = 2or3. 
It is considered a nonperturbative quantity by excellence, in the sense that it is trivially vanishing at any finite order 
of (ordinary) perturbative QCD in the chiral (massless quark) limit, as we recall in more details below. 

There is a long history of its determination from various models and analytic methods, the best known being the 
Gell-Mann-0 akes-Renner (GMOR) relation [1], relating the quark condensate to the pion mass, the decay constant 
F n , and the light current quark masses m Ut d, typically for the degenerate two-flavor case: 

F 2 ml = —(m„ + m d )(uu) + 0(m 2 q ) . (1.1) 

Nowadays the light current quark masses can be precisely extracted from lattice simulations [2] or spectral sum 
rules (see e.g. [3]), giving from using the GMOR relation above an indirect precise determination of the conden¬ 
sate. However, as indicated relation (1.1) is valid upon neglecting possible higher order terms 0{rn 2 q ). Indeed, the 
GMOR relation entails explicit chiral symmetry breaking from current quark masses. Since the condensate is more 
a property of the QCD vacuum in the strict chiral limit, it is highly desirable to obtain a possible ‘first principle’ 
determination of this dynamical quantity in the strict chiral limit, to disentangle from quark current mass effects. 
An early analytic determination was in the framework of the Narnbu- Jona-Lasinio (NJL) model [4] and its various 
extensions as a low-energy effective model of QCD, valid at least for the gross features of chiral symmetry breaking 
properties. In the NJL model, the condensate is evaluated in the strict chiral limit, or taking into account explicit 
breaking from small current quark masses, in the leading (large- N) approximation [5] as a function of the physical 
cutoff and other parameters of the model to be fitted from data. There are also related other attempts based on 
analytic methods, like Schwinger-Dyson equations [6, 7] typically. Phenomenological values of the condensate can 
also be extracted [3, 8] indirectly from data using spectral QCD sum rule methods [9], where the quark condensate 
and other higher-dimensional condensates enter as nonperturbative parameters of the operator product expansion in 
inverse powers of momenta. 

More recently ‘ab initio’ lattice calculations have determined the quark condensate by several independent ap¬ 
proaches, some actually related to the GMOR relation, or some more direct determinations with different methods [10] 
(see [2] for a review of various recent lattice determinations). In particular after early pioneering work [11], there 
has been more recently a renewed intense interest in computing on the lattice the spectral density of the Dirac 
operator [12, 13], directly related to the quark condensate through the Banks-Casher relation [14], However, while 
many recent lattice results are statistically very precise, lattice determinations rely in the end on extrapolations to the 
chiral limit, often using chiral perturbation theory [15] input for that purpose. Earlier general results on the spectral 
density in the nonperturbative low eigenvalue range had been obtained in [16], and attempts to use chiral perturbation 
theory information was elaborated e.g. in [17]. The link between different definitions of the quark condensate, in 
particular between the spectral density and the condensate appearing in the operator product expansion (OPE), has 
been carefully discussed in [7]. 

On more phenomenological grounds there are long-standing questions also on the dependence of the quark conden¬ 
sate on the number of flavors (n/ = 2 or n/ = 3 for the physically relevant cases). In particular it had been advocated 
that the GMOR relation may receive substantial corrections, and some authors indeed found significant suppression 
of the three-flavor case with respect to the two-flavor case [18], which may be attributed to the relatively large explicit 
chiral breaking from the not-so-small mass of the strange quark, roughly of order m s < A.qcd/3- There are also 
some hints that chiral perturbation theory might not converge very well for n/ = 3 [18, 19], for which value lattice 
simulations also show large discrepancies between different collaborations and methods (with some results with, and 
some other results without relative suppression of the n/ = 3 condensate [2]). 

Our recently developed renormalization group optimized perturbation (RGOPT) method [20-22] appears particu¬ 
larly adapted to estimate this quantity, since it gives a nonperturbative sequence of approximations starting from a 
purely perturbative expression. This allows by construction an analytic ‘first principle’ calculation of this quantity 
and gives a nontrivial result by construction in the chiral limit, in contrast with the standard perturbative one (see 
also [23, 24] for earlier attempts in that direction). Morever this also gives us a very simple analytic handle on the 
exploration of the chiral limit for arbitrary number of flavors (2 or 3 in practice), which in our framework is simply 
contained in the known flavor dependence of the first few perturbative coefficients, while this appears more difficult 
at present both for chiral perturbation theory and lattice simulations. 
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The paper is organized as follows. In section II we shortly recall basics of the spectral density and its Banks-Casher 
connection with the condensate. In section III we recall the main OPT method and our RGOPT version incorporating 
consistent RG properties. We adapt the RGOPT to the spectral density case in section III.C. Section IV is a digression 
where we first consider the spectral density RGOPT calculation in the Gross-Neveu model, where it can be compared 
with the exact result for the fermion condensate in the large -N limit, known from standard methods. Section V deals 
with the actual computation of the optimized spectral density in QCD at the three presently available orders (two, 
three and four loops) of the variationally modified perturbation. Detailed numerical results are presented as well as 
some comparison with other recent determinations of the quark condensate. Finally section VI is a conclusion. 


II. SPECTRAL DENSITY AND THE QUARK CONDENSATE 

We shall just recall in this section some rather well-known features of the spectral density and its connection with 
the chiral condensate, known as the Banks-Casher relation [14] (see also e.g. [16] for more details), to be exploited 
below. We thus start from the (Euclidean) Dirac operator which formally has eigenvalues X n and eigenvectors u n , 

i Ip u n (x) = X n u n (x); lp = 0 + g^, (2.1) 

where Ip is the covariant derivative operator and A the gluon field. Except for zero modes, the eigenvectors come 
in pairs {u n (x); 75 u n (x)}, with respective eigenvalues {A n ; —A ra } that depend on A. In the discrete case (be. on a 
lattice with finite volume V), by definition the spectral density is given by 

p(A)4(E^(A-AL- ] )), (2.2) 

n 

where S(x) here is the Dirac distribution and (• • •) designates averaging over the gauge field configurations, 

N 

() = / [dA] det(L0 + m). (2.3) 


The quark condensate is given by 


V Jv d4x ^ x ^ x ^ = ~ 2 y 


A„>0 


AS 


i2' 


Now when V goes to infinity the operator spectrum becomes dense, so that 

pW 


pOO 

{qq} = —2m dA 

J 0 


A 2 


m 


2 ’ 


(2.4) 


(2.5) 


where p{ A) is the spectral density. 

The Banks-Casher relation is the m —> 0 limit of this, giving the condensate in the relevant chiral symmetric limit 
as 


lim (qq) = -7iy>(0), (2.6) 

ra—>-0 

if the spectral density at the origin can be known. This is an intrinsically nonperturbative quantity, vanishing to all 
orders of ordinary perturbation, just as the left-hand side of this last equation. Now taking into account that for 
non-zero fermion masses m, (qq) = ( qq)(m ) = —E(m) we have from the defining relations (2.2), (2.5) the following 
interesting tautology, 

p( A) = 2- [E(iA + e) - E(iA - e)[ | e ^ 0 , (2.7) 

be. p( A) is determined by the discontinuities of E(m) across the imaginary axis. The interest of this relation is 
that, when the quark mass is nonzero, (qq) has a purely perturbative series expansion, known to three-loop order 
at present, and its discontinuities are simply given by those coming from the perturbative, purely logarithmic mass 
dependence. Therefore it makes sense to calculate a perturbative spectral density using the above relation. Usually it 
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will be of little use, since taking its A —> 0 limit, as relevant for the true chiral condensate, will only lead to a trivially 
vanishing result [7]. But the OPT (and in particular our RGOPT version) series modification after performing the 
variational 6 expansion (see below), is precisely the analytic handle giving a nontrivial result for A —> 0, just as it 
gives a nontrivial result for m —> 0 for modified perturbative series with mass dependence. 

Thus the recipe we shall apply is clear: In a first stage we calculate the purely perturbative expression of p( A) up to 
four loops, using the logarithmic discontinuities involved in Eq.(2.7). Then we perform a variational transformation 
(so-called 5-expansion, defined in next section) on the perturbative series, and solve appropriate optimization OPT 
and renormalization group (RG) equations (to be precisely defined in next section) to derive a nontrivial, optimized 
value of the relation (2.6). 


III. OPTIMIZED AND RG OPTIMIZED PERTURBATION (RGOPT) 

A. Standard OPT 

The key feature of the optimized perturbation (OPT) method (appearing in the literature under many names and 
variations [25]) is to introduce an extra parameter 0 < S < 1, interpolating between Cf ree and C lnt for any Lagrangian, 
in such a way that the mass parameter in is traded for an arbitrary trial parameter. This is perturbatively equivalent 
to taking any standard perturbative expansions in the coupling g(p), after renormalization in some given scheme ( e.g. 
MS -scheme with arbitrary scale m w), reexpanded in powers of 6 after substituting, 

m —> m (1 — 5)°, g —> 5 g . (3.1) 

Such a procedure is consistent with renormalizability [23, 24, 26] and gauge invariance [24], whenever the latter is 
relevant, provided of course that the above redefinition of the coupling is performed consistently for all interaction 
terms and counterterms appropriate for renormalizability and gauge invariance, as is the case for QCD. In (3.1) we 
introduced an extra parameter a, to reflect a certain freedom in the interpolation form, which will be crucial to 
impose compelling RG constraints, as discussed below and in our previous work [21, 22]. Applying (3.1) to a given 
perturbative expansion for a physical quantity P(m, A), reexpanded in S at order k, and taking afterwards the 6 —> 1 
limit to recover the original massless theory, leaves a remnant m dependence at any finite 5 k order. The arbitrary 
mass parameter m is then most conveniently fixed by an optimization (OPT) prescription, 

jLp( k \m,g,S = l)| m=A = 0, (3.2) 

am 

which generally determines a nontrivial optimized mass m(g), having a nonperturbative g dependence, realizing 
dimensional transmutation. (More precisely, for asymptotically free theories, the optimized mass is automatically 
of the order of the basic scale A ~ /j,e~ 1 ^ b ° 9 \ in contrast with the original vanishing mass). In simpler (D = 1) 
models the procedure may be seen as a particular case of “order-dependent mapping” [27], which has been proven[28] 
to converge exponentially fast for the D = 1 <f> 4 oscillator energy levels. For higher dimensional D > 1 renormalizable 
models, the behavior at large orders in 5 is more involved, and no rigorous convergence proof exists, although OPT was 
shown to partially damp the factorially divergent (infrared renormalons) perturbative behavior at large orders [29]. 
Nevertheless the OPT can give rather successful approximations for nonperturbative quantities beyond mean-field 
approximations in a large variety of models [25, 30, 31], including studies of phase transitions at finite temperatures 
and densities [33, 34]. 


B. Renormalization group optimized perturbation (RGOPT) 

In most previous OPT applications [25], the linear 6 expansion is used, namely assuming a = 1 in Eq. (3.1) mainly 
for simplicity and economy of parameters. However a well-known drawback of this conventional OPT approach is that, 
beyond lowest order, Eq. (3.2) generally gives more and more solutions at increasing orders, many being complex¬ 
valued, as a result of exactly solving algebraic equations in g and/or m. This problem is typically encountered first at 
two-loop order. In general, without some insight on the nonperturbative behavior of the solutions, it can be difficult 
to select the right one, and unphysical nonreal solutions at higher orders are embarrassing. As it turns out, RG 
consistency considerations provide a compelling way out, as developed in our more recent approach [20-22], which 
differs crucially from the more conventional OPT based on the linear 6 expansion in two main respects. First, it 
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introduces a straightforward combination of OPT and RG properties, by requiring the (5-modified) expansion to 
satisfy, in addition to the OPT Eq. (3.2), a standard RG equation, 


(^ (fc) (m, g,5 = 1)) = 0 , 


where the RG operator is defined as usual 1 , 


cl d d d 


(3.3) 


(3.4) 


Note, once combined with Eq. (3.2), the RG equation takes a reduced simple form, corresponding to a massless theory, 

P {k \m,g,6 = 1) = 0. (3.5) 

Thus Eqs. (3.5) and (3.2) together completely fix optimized m = m and g = g values. 


d o^ d 

^ + m T s 


Remark indeed that RG invariance is in general spoiled after the rather drastic modification from (3.1), reshuffling 
interaction and free terms from the original perturbative expansion. This feature has seldom been considered and 
appreciated in previous applications of the OPT based on the linear (5-expansion method to renormalizable theories. 
Thus RG invariance has to be restored in some manner, accordingly Eq. (3.5) gives an additional non-trivial constraint. 
Intuitively, just as the stationary point OPT solutions from Eq. (3.2) are expected to give sensible approximations, 
at successive orders, to the actually massless theory, one similarly expects that combining the OPT with the RG 
solutions of (3.5) should further give a sensible sequence of best approximations to the exactly scale invariant all 
order result. (N.B.: An earlier way of reconciling the S expansion with RG properties was used in [23, 24, 26]: 
Schematically it amounted to resumming the 5-expansion to all orders, which can be done in practice only for the 
pure RG dependence up to two-loop order. These resummations came as rather complicated integral representations, 
rendering difficult generalizations to higher orders, other physical quantities or other models of interest. In contrast 
the purely perturbative procedure of imposing Eqs. (3.2), (3.5) is a considerable shortcut, straightforward to apply 
to any model, being based solely on purely perturbative expansions.) 

Yet applying Eq. (3.2), (3.5) without further insight still gives multiple solutions at increasing orders. So we 
proposed [21, 22] a compelling selection criterion, by retaining only the branch solution(s) g(m) (or equivalently 
m(g)) continuously matching the standard perturbative (asymptotically free (AF) RG behavior in the QCD case) for 
vanishing coupling, namely, 

g(n » fh) ~ (26 0 In -^) _1 + 0 ((In -^-) -2 ). (3.6) 

m m 


Now the crucial observation is that requiring at least one of the solutions of Eq. (3.5) to satisfy (3.6) implies a 
strong necessary condition on the basic interpolation (3.1), fixing the exponent a uniquely in terms of the universal 
(scheme-independent) first order RG coefficients [21, 22]: 


_ 7o 
2&o 


(3.7) 


which is the second important difference of the present RGOPT with respect to the standard OPT 2 . For the critical 
value (3.7), Eq. (3.5) is in fact exactly satisfied at the lowest 5° order, therefore giving no further constraint. At 
higher 6 orders, (3.7) implies that one at least of both the RG and OPT solutions fulfills Eq. (3.6), and solutions with 
this behavior are essentially unique (although not necessarily) at a given perturbative order. Moreover, taking (3.7) 
drastically improves the convergence of the method: More precisely the known nonperturbative result of generic pure 
RG-resunmied expressions are obtained exactly from RGOPT at the very first 5-order, while the convergence of the 
conventional OPT with a = 1 is not clear or very slow, if any (see Sec. III.C of ref. [22] for details). 

The criterion (3.6) can easily be generalized to any model, even nonasymptotically free ones, by similarly selecting 
those optimized solutions that simply match the standard perturbative behavior for small coupling values. Thus 
clearly the resulting unique critical value like in (3.7) is valid for any model with its appropriate RG coefficients. For 


1 For QCD our normalization is /3(g) = dgjd 1 n // = —2bog 2 — 2b\ rf’’ + ■ ■ •, Zrn(f'j) = 70S + 7 1 .7“ + • ■ •, where g = Arras- The bi, 7 i 
coefficients up to four loops are given in [35]. 

2 The important role of the anomalous dimension 70 / (2&o) appeared also in our earlier constructions resumming the RG dependence of 
the (5-expansion [23, 24, 26, 29], although it had not been recognized at that time as a crucial RG consistency property of lowest 8 orders. 
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the QCD spectral density, as we will see below the equivalent of the criteria (3.6) indeed selects a unique solution at 
a given order for both the RG and OPT equations, at least up to the four-loop order (as was also the case for the 
pion decay constant [22]). 

Incidentally, a connection of the exponent a with RG anomalous dimensions/critical exponents had also been 
established previously in a different context, in the D = 3 $ 4 model for the Bose-Einstein condensate (BEC) critical 
temperature shift, by two independent approaches [30, 31], where for this model it also leads to real OPT solutions [31]. 
Indeed in ref. [30, 32] it was convincingly argued, based on critical behavior considerations, that the OPT can only 
converge if an appropriate Wegner critical exponent is used in the interpolation (3.1), which appears quite similar to 
our criterion (3.7). Note however that (3.7) is identified exactly from the known first-order RG coefficients, thus valid 
for any model, w r hilc in [32] the analoguous exponent was determined more approximately by looking for a plateau in 
the variational parameter dependence. In any case, from these examples, it is established that it is necessary for the 
OPT method to give useful results to have in general an exponent a in (3.1) differing from 1 in a well-defined way. 

Coming back to the present OPT and RG Eqs. (3.2), (3.5), beyond lowest orders AF-compatible solutions with 
behavior (3.6) are however not necessarily real in general. A rather simple way out is to further exploit the RG 
freedom, considering a perturbative renormalization scheme change to attempt to recover RGOPT solutions both 
AF-compatible and real [22]. In the present case of the spectral density, this extra complication is not even necessary, 
at least up to the highest order here studied (four loops): We shall find also that the unique AF-compatible RG and 
OPT solution remains real. 


C. RGOPT for the spectral density 


Formally a generic pertubative expansion for the condensate reads typically 

(qq) pert = ™ 3 9 P H fp k lnP_fe (~) > ( 3 - 8 ) 

p k =o ^ 


where the f p k coefficients are determined by RG properties from lowest orders for k < p. According to Eq. (2.7), 
calculating the (perturbative) spectral density formally involves calculating all logarithmic discontinuities. This is 
conveniently given by taking in any typical perturbative expansion like (3.8), all nonlogarithmic terms to zero, those 
trivially having no discontinuities, while replacing all powers of logarithms, using m —> i|A| etc., as 


In" 


= — In" 
2 n 


m 


P 


1 1 

2" 2br 


2 In 


JA[ 

P 


leading to the following simple substitution rules for the first few terms 

|A| 


In 


m 




1/2; In 


In ■ 


In' 3 


P 


— ( 2 In -—- — i7T 
P 


- In 2 — 

2 /i 


(3.9) 


(3.10) 


and so on (note the appearance of nonlogarithmic ~ tt 2 terms starting at order In 3 to). This gives a perturbative 
expression of the spectral density of the generic form 


Ppert (\X\,g) = |A| 3 ^ g p ^2 fpk ln p_fc (—), 

p>l fc=o ^ 


(3.11) 


where the determination of the coefficients follows from the above relations (3.9). 

To obtain the RG equation for p(g, A), we use the defining integral representation of the spectral density in (2.5) and 
the basic algebraic identity 


dm. d A 

dm A 2 + m 2 dX A 2 + m 2 ’ 


(3.12) 


Throwing away surface terms in partial integrations as usual in the spirit of dimensional regularization, one thus finds 
that p{ A) actually obeys the same RG equation as {qq ), with dm. replaced by dX, 




p(g, A) = o. 


(3.13) 
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One can next proceed to the modification of the resulting perturbative series p( A, g) as implied by the (5-expansion, 
now, from Eq. (3.12) clearly applied not on the original mass but on the spectral value 3 A: 

A -+ A(1 — S) a g^6g. (3.14) 


Optimizing perturbation theory means that the derivative with respect to to of 


OO 


E 

n—0 


(-i) r 



(3.15) 


being formally zero 4 , one should obtain a good approximation for the value at to = 0 of ( qq ) at finite order by setting 
to zero the derivative of a finite number of terms of this series, see Eq. (3.2). Using Eq. (3.12), this mass optimization 
on (qq) thus translates into an optimization of the spectral density with respect to A, 

- ^=0, (3.16) 

at successive S k order. 



IV. LESSONS FROM THE GROSS-NEVEU MODEL 


The fermion condensate can also be defined from the spectral density for the .D = l + 1 O(N) Gross-Neveu (GN) 
model [36]. This will give us a very useful guidance for the more elaborate QCD case below, and we also set up some 
formulas actually generically valid for both the GN model and QCD. 


We start from the known expression of the vacuum energy evaluated in the large -N limit [26] after all necessary 
mass, coupling, and vacuum energy (additive) renormalizations, in terms of the explicit mass m and mass gap M(m, g), 


( N 

Egn = " t 

where the mass gap is defined in compact form as 


to . 


M 2 (m, g) + 2 —M (to, g) 


M 

M (to, g) = to I 1 + g In — 


-l 


(4.1) 


(4.2) 


where to = m(/z) and g = g(p) are the renormalized mass and coupling in the MS scheme (after convenienly rescaling 
the original coupling defined by (l/2)gg. JV (4 , 4') 2 as gQ N N/ir = g). The fermion condensate is formally given by the 
derivative with respect to to of the vacuum energy, giving simply after some algebra 

V ipip)GN(m,g ) = - I — I --- • (4.3) 

This means that up to a trivial overall factor, the fermion condensate is directly related to the mass gap, as intuitively 
expected also from mean-field arguments. Eq. (4.3) has a well-defined nontrivial perturbative expansion to arbitrary 
order, 

^ (^)gjv(to, g) = — (5E , \E , )Giv — to — L m + gL m (L m + 1) — —L m (L m + 2)(2 L m + 1 )g 2 + 0(g i )^ j(4.4) 

where L m = him/p (and for convenience we redefined the condensate by a trivial rescaling). From properties of the 
implicit M(m) defined by (4.2) and its reciprocal function m(M), one can establish [26] that M(m) —> A = fie _1 / 9 
for to —> 0 , which translates here into the simple relation 

-(#<J>)o/v(to -+ 0) = — , (4.5) 

9 



3 We adopt in the following the notation A = |A| since it is necessarily positive. 

4 For simplicity, we have set a to one in this equation. 



which provides a consistent bridge between the massive and massless case. But deriving this needs the knowledge of 
the all-order expression (4.2), only known exactly in the large- N limit. 

Now alternatively, performing the substitution 5 (3.1), expanding at order 5 P , setting 6 to one, and optimizing the 
resulting expression with Eqs. (3.2) and (3.5) gives the exact result (4.5) at any order in S p , at optimized coupling 
and mass values, 

9= 1; L m = - 1, (4.6) 

just as in the mass gap case [20]. This is not very surprising, in view of the rather trivial relation (4.5) between the 
mass gap and the condensate in the large- N limit. 

However here we shall try to obtain this large-TV result in an indirect way, using the spectral density, the aim being 
evidently to test on the exactly known result (4.5) the possibility of calculating a spectral density and of estimating 
its reliability from the first few perturbative orders, in order to subsequently apply the same procedure in the more 
challenging QCD case. 

Thus from (4.4) and using m -+ i|A| and relations (3.9) the perturbative expression of the spectral density, for 
instance restricted up to four-loop order g 3 , follows as 

Pgn (^> 9) — 2 + 9 + 2 ^ — + 20La + 4 — n 2 ) + 

(48 L\ + 156L| + (108 - 12t t 2 )L x + 12- 13tt 2 ) + 0(g 4 )^j , (4.7) 

where now L\ = In X/p. We can then proceed by applying (3.14), expanded to order S k , then taking S —> 1, and finally 
applying on the resulting expression the OPT (3.16) and RG (3.13) Eqs. In practice we shall proceed to relatively 
high perturbative orders, since the exact expression (4.2) may be formally expanded to arbitrary order. 

Thus with those obvious replacements we can proceed first with the (5-expansion (3.14) operating on A and g , and 
taking the relevant value of the exponent for the large- N case, a = jo/(2bo) = 1 in (3.1). At first nontrivial (5 1 order 
we simply obtain 


p s (6 —> 1; A, g) 



from which the OPT (3.16) and RG (3.13) give the unique solution, 


Lx 





(4.8) 


(4.9) 


which plugged back into (4.8) gives the final result, using also the Banks-Casher relation (2.6), 

9 {^^)Gw{rn -+0)/A = -7rgp(A, g)/A = ~ yfe ~ -1.2949 , (4.10) 

which is to be compared to the exact result g (T 4 , ) ( 3 jv(to —> 0)/A = —1 in this normalization. We then proceed to 
rather high perturbative order, which in practice becomes relatively involved algebraically but can be easily handled 
with computing softwares like Mathematica® [37]). 

The results up to order g 18 (beyond which numerics become really tedious) are given in Table I (where we give for 
convenience the scale invariant condensate 5 (TT) values). Actually, starting at order g , there are more than one 
real solution. We show here those solutions which are unambiguously determined to be the closest to the standard 
AF perturbative behavior, L\ ~ —l/g + 0(l) for g —> 0, by analogy with the exact value —gL m = 1 obtained for the 
simpler mass gap case. One can see a regular pattern for the values of the optimized coupling for such solutions, with 
g ~ 0.4, while the optimal spectral parameter L\ is less stable with wider variations at successive orders. 

These approximate results at successive orders when optimizing the spectral transform are clearly in contrast 
with those obtained from optimizing directly the original perturbative expansion, where as explained above the AF 


Taking a = 1 as appropriate for the large-AT GN model. 


5 



9 


TABLE I. RGOPT Si^loN an( j corresponding optimized coupling g and (logarithm of) optimized eigenvalue lnA/p results 
at successive orders in S. 


S k order 


9 

In — 

M 

A 

1 

1.295 

1 

2 

3 

2 

2 

0.984 

0.398 

-2.547 

3 

0.897 

0.335 

-3.278 

4 

1.081 

0.399 

-1.373 

5 

0.924 

0.394 

-1.919 

6 

0.877 

0.372 

-2.337 

7 

0.980 

0.386 

-1.332 

8 

0.903 

0.389 

-1.716 

9 

1.065 

0.358 

-0.928 

10 

0.934 

0.383 

-1.312 

11 

0.894 

0.385 

-1.613 

12 

0.978 

0.364 

-1.001 

13 

0.972 

0.349 

-3.77 

14 

1.013 

0.339 

-4.049 

15 

0.989 

0.403 

-2.902 

16 

1.027 

0.391 

-3.107 

17 

0.978 

0.428 

-2.405 

18 

1.013 

0.420 

-2.574 


compatible solution always gives the exact result for the condensate, with the corresponding optimal coupling and 
mass values (4.6) obtained already at first and all successive orders. Here we see in Table I that the second order 
is very close, less than 2%, from the exact result, but at higher orders the solutions exhibit a rather slow empirical 
convergence, with an oscillating behavior towards the exact large -N limit result. 

This slow convergence can be essentially traced to the effect of numerous factors ^ r 2p , p = 1, ...n/2 appearing at 
perturbative orders g n , g n+1 for n > 2 from the discontinuities. These terms clearly spoil the originally neat simple 
form of the large-TV resunnned mass in (4.2). Moreover, starting at order g 3 these n 2p terms come with larger 
coefficients relative to the other non-logarithmic coefficients, originating from the In rn coefficients in the original 
perturbation (4.2), (4.7) which are roughly all of the same order 0( 1). More precisely inspecting (4.7) one can see 
that at order g 2 the —tt 2 /8 contribution is numerically very roughly twice as large as the other non-logarithmic 
coefficient (1/2), while at order g 3 the relevant terms to be compared are the two last ones: 12/24 and —(13/24)7r 2 , so 
the latter is an order of magnitude larger than the former. This explains the very good result at order g 2 and also the 
degraded result at the next g 3 order. A similar behavior is observed at higher orders, until at sufficiently high order 
the original perturbative coefficients of In to also start to grow fastly, such that a balance with the n 2p contributions 
can again occur, with the (slow) convergence observed. Indeed the spectral parameter optimization (3.16) tends to 
damp these relative n 2p contributions: For instance the combined contributions of all L>-dependent terms for the 
optimal value L\ ~ —3.28 at order g 3 in Table I almost cancel the large 7r 2 term. But since the latter terms are 
growing with the order, it is not surprising that the optimal L\ values are not much stable in Table I. 

Interestingly however, if instead of taking the exact results at a given order, we consider a well-defined approxima¬ 
tion, by keeping, at growing 6 n order, only tt 2p terms with a fixed maximal power, then there is a higher S n order at 
which one recovers the simple exact RGOPT solution 6 : g = 1,L\ = —1, ('k'P )gn = —A. For instance, keeping only 
7T 2 and 7r 4 terms (the latter appearing first at order g 4 ) and increasing the S k order, the exact solution is recovered at 
order 6 s . This cancellation mechanism thus indicates that the “maximal convergence” properties of RGOPT, specific 
of the simpler GN model mass gap expression in (4.2) [20], are not completely lost within the perturbative spectral 
density, but rather hidden, being obstructed by the more involved perturbative coefficients. Those remarks are to be 
kept in mind when comparing with the QCD case below. 


This nice property may appear somewhat artificial but it can be explained more rigorously: As observed in [20] the mass gap M(m,g) 
at a given (sufficiently high) order S n , has flat optima roughly at order ~ n/2 ( i.e. its n/2th derivative with respect to m vanishes). 
Now for the spectral density, 7r 2p terms appearing at order g p arise from (2.7) as the coefficient of the (logarithmic) derivatives d/d In m 
of order p : So if discarding terms of higher power n 2q ,q>p , there is necessarily a fixed higher order at which the cancellation of all 
7 r 2p terms occur. 
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V. DETERMINATION OF THE QCD QUARK CONDENSATE 
A. Perturbative three-loop quark condensate 



FIG. 1. Samples of standard perturbative QCD contributions to the chiral condensate up to three-loop order. The cross denotes 
a mass insertion. 


The perturbative expansion of the QCD quark condensate for a nonzero quark mass can be calculated systematically 
from the related vacuum energy graphs. A few representative Feynman graphs contributions at successive orders up 
to three-loop order are illustrated in Fig. 1 (there are evidently a few more three-loop contributions not shown here). 
Note that the one-loop order is 0(1) = O(g 0 ). The two-loop contributions were computed long ago [38] and the 
three-loop ones in [40]. The three-loop order result in the MS'-scherne reads explicitly 


m ( qq)QCD (m, j) = ^m 4 ( ^ - L 


_l_ JL(t,2 _ ^ 

i o \^m a -kn 


7T 


5 g 
+ 12 ) + ( 16^ 


) 2 93 (A ,n f ) 


(5.1) 


where m = and g = 47 rag(/Li) are the running mass and coupling in the MS scheme, and the three-loop coefficient 
reads 7 [40] 

< 73 (A, n f ) = y (6185 - 768 a 4 - 32 In 4 2 + 192 In 2 2 z 2 + 504 2 3 + (672 2 3 - 750) n/ + 528 z 4 ) 

+ ( 52 n f - yy + y ^ 3 )Lm - y (5 n f - 141)L^ + y (2 n f - 81 )L 3 m , (5.2) 

where Zi = £(i) and a 4 = Li 4 (1/2) . 

The calculation in dimensional regularization of (5.1) actually still contains divergent terms needing extra subtrac¬ 
tion after mass and coupling renormalizations in the MS scheme. The correct procedure to obtain a RG invariant 
finite expression when subtracting those divergences consistently is well known in the standard renormalization of 
composite operators with mixing [39]. We can define [24] the needed subtraction as a perturbative series, 


sub ( 5 , m) = —^2 

n ^' 


Sig 


i> 0 


(5.3) 


with coefficients determined order by order by 

= Remnant( 0 , m) = fi-^-[m(qq)(pert) \ finite], (5.4) 

d/i d/i 

where the remnant part is obtained by applying the RG operator Eq. (3.4) to the finite expression (5.1), not separately 
RG invariant. Eq. (5.3) does not contain any lnm//i terms and necessarily starts with a so/g term to be consistent 
with RG invariance properties. To obtain RG invariance at order g k , fixing Sk in (5.3), one needs knowledge of the 
coefficient of the lnm term (equivalently the coefficient of 1/e in dimensional regularization) at order g k+1 . Concerning 
the condensate the extra contribution to the RG equation (5.4) is the so-called anomalous dimension of the QCD 
(quark) vacuum energy, entering the renormalization procedure of the m(qq) operator due to mixing with m 4 x 1, 


' The originally calculated expressions in refs. [40] are given for arbitrary N c colors, massive quarks and n; massless quarks entering 
at three-loop order. In our context m is the (nj -degenerate) light quark mass and its precise mass dependence is what is relevant for 
the optimization procedure, so one should trace properly the full nj, n/ t dependence, and take n; = 0. = ri j with nf = 2 (3) for the 

SU(2) (resp. SU( 3)) case. 
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and also given explicitly to three-loop order in [40]. The s* coefficients can be expressed in terms of RG coefficients 
and other terms using RG properties. In compact form for completeness in our normalization they read 8 , 


so 


si 


S2 


S3 


1 

47t 2 (6 0 - 270) ! 

1 1 bi — 2 71 

12-7T 2 4 bo — 2 70 

(112077 + 24519 n f + 2101 n) + 576(15 - 58 n f ) z 3 ) 

1152 7r 4 (—81 + 2n/)(15 + 2 n f ) ’ 

4.710 _lo a 4 (—81 + 2n/)(—4.615 + n/)(15 + 2 n f ) + 2.47 + 1.275n/ + 0.303n^ 

(-57 + 2n/)(—81 + 2n/)(15 + 2n/) 


0.019n^ + 3.110- 4 n) 


(5.5) 


B. Perturbative spectral density 


One crucial advantage of using the spectral density with the Banks-Casher relation (2.6) is that it gives a direct 
access to the QCD condensate in the chiral limit, unlike the original direct RG invariant expression m(qq) in (5.1), 
where the condensate is being screened by the mass for m —> 0. Indeed, note that a direct RGOPT optimisation in 
the GN model of the corresponding expression m(4 , 4/) gives an exactly vanishing result consistently at any orders [22] 
(even though the optimized GN mass is clearly nonvanishing, rh = A). 

Taking thus the logarithmic discontinuities according to Eqs. (2.7),(3.9) gives us the perturbative spectral density 
up to three-loop order, 


PQCdW ~ 2^2 (“S + *2 ( iA 



(5.6) 


where now L\ = ln(|A|//i) and 

4 D {\rif) (52n/ - + y z 3 ) - y(5 n f - 141 )L X 

+|(2n / -81)(^|-^). (5.7) 

Note that the 7 r 2 in the last term arises from the discontinuities of ln 3 (m 2 ) according to (3.9). Now remark that none 
of the nonlogarithmic contributions in the original perturbative expression (5.1) contribute to the spectral density. 
Thus similarly all subtraction terms in (5.5), which are necessary for RG invariance of the original expression, do not 
contribute as well, thus making the final expression to be optimized relatively simpler. 

This point is worth elaborating in some detail. Just as for the GN model, instead of using the spectral density we 
could apply the RGOPT method more directly to the original perturbative expression of the condensate, Eq. (5.1), 
including in this case the subtractions (5.3), (5.5) required by RG invariance (and also removing an overall factor m 
from Eq. (5.1) to define a nontrivial (qq) in the chiral limit). When this is done, one finds rather unstable values for 
the optimized mass, coupling, and resulting condensate, showing no clear empirical convergence pattern at successive 
orders, at least at the presently available (three-loop) order. Furthermore these results tend to give a wrong sign 
(positive, or ambiguous) condensate. More precisely, at the first nontrivial 5° (one loop) order, firstly there is no 
common nontrivial RG and OPT Eqs. solution. Considering then the OPT or RG Eqs. alone, both give a positive 
condensate, of roughly the right order of magnitude: (qq)(S°) = y/e/{2TT 2 )h.\ rs ~ 0.08A|— from the OPT, and a very 
similar value is obtained from the RG. Next at the (5 1 (two-loop) order, the (unique) AF-compatible branch solution 
of the combined RG and OPT Eqs. (3.5), (3.2) gives complex-valued optimized coupling, mass and condensate, 
with a negative real part for the condensate but a much larger imaginary part: (gg)(5 1 ) ~ (—0.08 ± 0.37 t)A|^-, 
a result clearly ambiguous. These calculations are also not very stable upon different truncations of perturbative 
higher order terms in the RG equation (3.5). Furthermore, attempting to recover real AF-compatible solutions by 
a perturbative renormalization scheme change, both at orders 6 and 6 2 , happens to give no solutions (unlike for the 


8 The expression of S 3 here approximated to 10 y ‘ relative uncertainty, largely sufficient for our purpose, needs the knowledge of the 
four-loop In m/ji coefficient, or alternatively the four-loop vacuum energy anomalous dimension. The latter, not explicitly available in 
the published literature so far, has been kindly provided to us by K. Chetyrkin and P. Maier [41] from a related work. 
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pion decay constant case where the imaginary parts were small enough to allow for such a scheme change with very 
stable results [22]). 

We can trace this wrong sign and unstable behavior to the fact that in four dimensions the (presumably dominant) 
one-loop contribution to the fermion condensate, given by the very first graph in Fig. 1, is quadratically divergent. 
The contribution of this quadratic divergence has actually the correct negative sign. Incidentally, in the Nambu- 
Jona-Lasinio model [4], an effective cut-off handles this divergence, and the quark condensate has automatically the 
correct sign. Note that in the NJL model the condensate (or equivalently the mass gap in the widely used leading 
order large -N approximation) is precisely given by the very same one-loop first graph of Fig. 1, up to trivial overall 
factors, while genuine QCD contributions only enter at the next orders with gluon and further quark loop dressing. 
In dimensional regularization, at lowest orders in the coupling, one finds that the extra subtractions (5.3), (5.5) have 
a sign opposite to the sign of the similar terms in the GN model, and that this is due to the fact that the pole 
of T(1 - D/2) in the perturbative calculation of the condensate changes sign when going from dimension D = 2 
(corresponding to the logarithmic divergence in the GN model, and quadratic divergence in the D = 4 NJL model) 
to dimension D = 4 (corresponding to QCD). Indeed, as is well-known most of the phenomenological successes of 
the NJL model rely strongly on the physical cutoff interpretation of the quadratically divergent mass gap in four 
dimensions 9 . 

Hence, the RGOPT appears bound to lead to a wrong sign and/or unstable results, if applied directly to the QCD 
perturbative expression of the quark condensate evaluated in dimensional regularization and related MS scheme at 
low orders. We do not know whether higher orders would cure this problem by ultimately stabilizing the result, 
but this appears rather unlikely since the first few orders are likely to remain dominant in our approach. (Indeed 
a general property of the optimized perturbation is that the optimized coupling g turns to be reasonably small, so 
that the first few orders dominate). We did not have this problem in our previous works [20-22] in which we were 
dealing with the pole mass and the pion decay constant, which are only logarithmically divergent quantities. Yet one 
should not hastily conclude that the OPT or RGOPT approaches are bound to fail in any situation where quadratic 
divergergences would be present in a cutoff regularization 10 . Rather, the above problems stress that in a given model 
it is crucial to choose carefully the basic entity to be perturbatively modified and optimized within the RGOPT 
framework. (This is analoguous to the traditional variationnal Rayleigh-Ritz method in quantum mechanics, where 
the trial wave functions should often be appropriately chosen to obtain a sensible result). This is why for QCD 
one must use instead the spectral density, which in our framework we anyway derived from the very same original 
perturbative condensate expression (5.1), but which at the same time formally gets rid of the influence of quadratic 
divergences. Indeed, only the infrared part A —> 0 in Eq. (2.5) can generate a nonzero result in the chiral limit, which 
is thus insensitive to ultraviolet divergences [16]. We note also that lattice evaluations of the condensate also bypass 
this potential quadratic divergence problem by using the spectral density [12, 13], or by extracting the condensate by 
more indirect methods, e.g. relying on the GMOR relation. 

We thus proceed with the actual RGOPT calculations for the spectral density at successive perturbative orders. 
First remark that, since there is no logarithmic L\ contribution in the spectral density at one-loop order (the one-loop 
In to contribution in (3.10) only giving the constant 1/2), there is no nontrivial A ^ 0 optimized solution of (3.16) at 
one-loop order. Thus we should start applying our method at the next two-loop order. 


C. Two-loop C>(<5) results 


Let us perform step by step the RGOPT optimization by restricting first (5.6) at the first nontrivial two-loop order. 
Concerning the J-expansion given by (3.14), it is crucial [21, 22] to take the right value of the exponent a, determined 
by the lowest order anomalous mass dimension, which makes the J-modihed series compatible with RG properties 
and matching asymptotic freedom (AF), as we recalled in some details in Sec. III.B above. In the case of the large -N 
limit of the GN model, one has simply a = yo/(26o) = 1. Actually, since m.(qq) is RG invariant to all orders rather 
than ( qq ), it is easily derived that the correct value to be used for {qq), and thus for the related spectral density from 
(2.5), is 


a = 


4,7o s 

3 { 2b 0 ' 


(5.8) 


9 The NJL model may be formulated in dimensional regularization to some extent, in dimensions 2 < D < 4, but with rather odd 
properties, see e.g. [42] 

10 In fact the (standard) OPT has been applied in the framework of the effective NJL model with a cut-off at next-to-leading J order [34], 
giving sensible results beyond the large-AT approximation, and consistent with important basic properties like the Goldstone theorem, 
the GMOR relation, etc. 
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Then to first nontrivial order in 5 the modified series reads, 


~Pqcd - ^2 A 




(5.9) 


and the OPT (3.16) and RG (3.13) equations have a unique solution, using also (2.6), given in the first line of Table II. 
For a simpler first illustration we actually used the RG Eq. (3.13) at the very first order with the one-loop coefficient 
bo, in order to get simple analytic solutions. Therefore we obtain ( qq ) 1 ^ 3 (n/ = 2)(/r ~ 2 . 2 A 2 ) ~ — 0.96 A 2 , a fairly 
decent value given this lowest nontrivial order. At two-loop perturbative order expression (5.9) does not depend 
explicitly of the number of flavors rif, but a rif dependence enters into the optimized results indirectly from the 
RG Eq. (3.13) involving bo(rif), also entering A„,. The corresponding optimized coupling is as = g/{ 47 r) ~ 0.83, a 
moderately large value very similar to the optimal coupling values obtained, at first nontrivial RGOPT order, when 
considering the pion decay constant F\ in ref. [22]. Of course the precise number obtained for the condensate depends 
on the precise definition of the A reference scale, which is generally perturbative and a matter of convention. To get 
the numbers in the first lines of Table II we have used the simpler one-loop form, A = /^e _1 ^ 2bo9 \ consistently with 
the one-loop RG Eq. used. When comparing below more precisely with other phenomenological determinations of the 
condensate, we will use a more precise perturbative definition of A, at four-loop order, in agreement with most other 
present determinations. Remark also that the condensate being scale dependent, our RGOPT optimization fixes also 
a scale, consistently with a defining convention for A, as indicated in Table II. 

For rif = 3 at order d one finds similarly the optimized results given in the first line of Table III, indeed very close to 
the rif — 2 results above. 

Now, since our basic expression originated from an exact two-loop calculation, it is a priori more sensible to apply 
the RG Eq. (3.13) at the same two-loop order, in order to capture as much as possible higher-order effects. Doing 
this we obtain the results given in the second lines of Tables II, III for rif = 2,3. Those results are therefore to 
be considered more consistent at two-loop order. One can already observe the substantial decrease of the optimal 
coupling as to a more perturbative value, and the correspondingly higher optimal scale /i, with respect to the results 
using pure one-loop RG equation. 

For completeness and later use we also give in Table II, III the corresponding values of the RG invariant condensate 
(qq)rgi , perturbatively defined in our normalization as 

(qq)rgi = (qq){v) (26 0 s)^° ( X + “ ^r) 9 + 0{g 2 )j , (5.10) 

where higher order terms not shown here are easily derived from integrating exp [J dg'y m (g)/(3(g)\ thus know pertur¬ 
batively to four-loop g 3 order since only depending on the RG function coefficients [35] bi, 7 *. The factor multiplying 
the scale-dependent condensate ( qq)(g ) in (5.10) is obviously the inverse of the one defining similarly a scale-invariant 
mass, given explicitly to four-loop order in the literature (see e.g. [43]). We will calculate the RG invariant condensate 
at successive orders in Tables II, III using (5.10) consistently at the same perturbative order as the RG order used 
in Eq. (3.13), and taking g = g = Anas, the corresponding optimized values obtained at each order. (Alternatively 
we could optimize directly the expression (5.10) instead of (qq)(p) (with the last term oc 7 m (g) in the RG Eq. (3.13) 
removed): This gives the same optimized solutions, as expected since it is formally completely equivalent, up to tiny 
numerical differences due to perturbative reexpansions, of less than 10 ~ 3 relative to the numbers given in Tables II, 
III.) 


TABLE II. Main optimized results at successive orders for n/ = 2, for the optimized spectral parameter A, the optimized 
coupling as, and resulting optimized condensate. We also give the RG invariant condensate (qq^ai calculated at the consistent 
perturbative order from (5.10). A 2 is conventionally normalized everywhere by Eq. (5.11), except in the very first line where 
the one-loop expression A = g e -1 /! 26 ° is rather used. 


5 k , RG order 

In — 

M 

as 


JL 

a 2 

-WtRGI 

A 2 

<5, RG 1-loop 
<5, RG 2-loop 

2275 

10092 

-0.45 

§f ~ 0.83 
0.480 

0.962 

0.822 

2.2 

2.8 

0.996 

0.821 

5 2 , RG 2-loop 
S 2 , RG 3-loop 

- 0.686 

-0.703 

0.483 

0.430 

0.792 

0.794 

2.797 

3.104 

0.792 

0.783 

<5 3 , RG 3-loop 
S 3 , RG 4-loop 

-0.838 

-0.820 

0.405 

0.391 

0.793 

0.796 

3.306 

3.446 

0.774 

0.773 
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TABLE III. Same as Table II for n/ = 3 


5 k order 

InA 

as 

-<-> 1/3 (a) 

A 

A 3 

-<59>R<77 

A 3 

S, RG 1-loop 

S, RG 2-loop 

283 

972 

-0.56 

m - °- 82 
0.474 

0.965 

0.799 

2.35 

3.06 

0.987 

0.789 

5 2 , RG 2-loop 
S 2 , RG 3-loop 

-0.766 

-0.788 

0.493 

0.444 

0.776 

0.780 

2.942 

3.273 

0.772 

0.766 

S 3 , RG 3-loop 
<5 3 , RG 4-loop 

-0.967 

-0.958 

0.414 

0.400 

0.769 

0.773 

3.540 

3.700 

0.745 

0.744 


D. Three-loop 0(5 2 ) results 


At three-loop g 2 order the n/ dependence enters explicitly within the perturbative expression of the spectral density, 
see Fig 1 and the last g 2 coefficient in Eq. (5.6), (5.7). This is interesting also in view of other results on the variation 
of the condensate value with the number of flavors [ 2 , 18]. 


We find a unique real AF-compatible optimized solution. More precisely, at this three-loop order there are actually 
two real optimized solutions for L\, a $, but the selection of the right physical solution is unambiguous since only one 
is clearly compatible with AF behavior for g = Anas —> 0, ln(A//x) ~ —dk/(2bog) + 0(1) with dk = 0(1), both for 
the RG and OPT equations . 11 In contrast the other real solution has for j->0a coefficient of opposite sign to AF, 
and gives La = lnA//x > 0, which also means incompatibility with perturbativity, since we expect /i A just like 
the perturbative range being g rh ^ A for the original expansion with mass dependence. Explicitly we obtain for 
n/ = 2,3 the results given in the third and fourth lines of Tables II, III respectively. More precisely as indicated there 
the third line results were obtained by taking the RG Eq. (3.13) truncated at two-loop order, while the fourth lines 
were obtained by taking the full three-loop RG equation. These results exhibit a very good stability when confronted 
with the relative arbitrariness in the order of the RG equation. Of course it is more legitimate to use the three-loop 
RG equation consistently at this three-loop order, that we shall adopt for our final determination of the condensate. 
Notice also that the optimal values as decreased by almost a factor two with respect to the lowest nontrivial order 
result above, which indicates that the resulting series is much more perturbative. But as almost does not change 
as compared to the more consistent two-loop results. Similarly as further slightly decreases and the optimal scale ji 
increases when going from two to three loops in the RG Eq. (3.13). 


As a matter of numerical detail, to obtain the results in Tables II, III we took a convention of the QCD scale A 
based on a perturbative four-loop expression [44]: 


A 4 — loop 

rif 


(g) = fie 


(bo g ) 


oh'2 

o exp 



b l ) + (R. 

b 2 ’ + [ 2b* 


bib-2 

b 2 o 



(5.11) 


This is convenient and important to compare precisely below with most recent other determinations, using the same 
four-loop perturbative order conventions for A. Actually when performing the RG Eq. (3.13) at order S k it would be 
more natural to adopt a A convention at the consistent (k + l)-loop order Afc +1 , given from (5.11) by taking 63 = 0 
(and 62 = 0) respectively at three loops (two loops). But this only affects an overall normalization of the final result, 
as A itself is not involved in the actual optimization process when performing Eq. (3.2) and (3.13). Besides, starting 
at three-loop order the differences obtained from such different conventions are minor. (The A convention also affects 
the precise value of the optimal scales fi /A in Tables II, III from which we shall start evolution to higher scale to 
compare with other determinations in the literature, see below). Strictly speaking, the different values of ( qq)(ii) 
obtained in Tables II, III e.g. at three-loop order cannot be directly compared, being obtained at different scales ft. 
Thus we also give the scale-invariant condensate values ( qq)RGi which can be more appropriately compared. 

Notice also that in spite of the more than 10% change in the optimal coupling as when taking two- or three-loop 
RG, the final physical value of the condensate only varied by 0.25% 12 : This also reflects a strong stability. Moreover 


11 Due to the nontrivial relations between lnm and InA, Eqs. (3.9), (3.10), the l/g coefficient of the correct AF branch L\(g) for g 0 
is not exactly — l/(2f>o), like it is for In rnf ft (Eq. (3.6)), essentially determined by the leading logarithms g k I)A(m///.) behavior. But 
this l/g coefficient has the correct AF sign, being at order S k : — <L,/(2&o) with dy, > 0 a. constant close to 1, which approaches slowly 
<ffc —> 1 as the perturbative order k increases. 

12 This result appears so stable partly due to the 1/3 power of the condensate in Tab. II, III. For the actual optimization performed on 
the quantity (qq ), the corresponding variation is ~ 0.7%. 
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the value of (qq ) 1 ^ 3 /A changed by about 20% with respect to the crude two-loop order result (first lines of Tables), 
but much less when compared to the more consistent two-loop order result (second lines). This shows a posteriori that 
stability appears at the first nontrivial two-loop result, with an already quite realistic value. This stability, suggesting 
that we remain within the domain of validity of perturbation theory, is an important requirement for the usefulness of 
our method. A similar behavior was observed when optimizing the pion decay constant in [22]. Remark also that the 
optimized coupling values as at successive orders happen to be rather close to those for which the scale invariance 
factor in (5.10) multiplying (qq)(p) would be exactly one (which for n/ = 2 happens for as — 0.483,0.461,0.457 
respectively at two-, three- and four-loop order). In other words, as is close to a (variational) “fixed-point” scale- 
invariant behavior. Had we found optimized as values a factor 2-3 smaller, or larger, we would obtain no valuable 
results beyond ordinary perturbation in the first case, or much more unstable results in the second case. However we 
stress that the optimal coupling as or optimal mass m do not have really a universal physical interpretation since 
the precise as and m values depend on the physical quantity being optimized. For instance when optimizing F n 
in ref. [ 22 ] at a given perturbative order, the corresponding as values were pretty much close to the present ones, 
but nevertheless slightly different. The physically meaningful result is obtained when replacing as and A within the 
quantity being optimized, like here the condensate. 

Comparing Tables II, III it also clearly appears that the ratio of the quark condensate to A 3 has a moderate 
dependence on the number of flavors n/, although there is a definite trend that (qq) rl ' :i is smaller by about 2 — 3% 

with respect to {qq)„ = 2 , in units of A n/ , at the same perturbative orders. The smallness of this difference was 
quite expected, due to the n/ dependence only appearing at three-loop order and the overall stability of the modified 
perturbation. However from various different estimations, including lattice [46], and ours [22]), there are some 
indications that A 2 > A 3 (although unclear from uncertainties, due to a larger uncertainty on A 2 ), which therefore 
could indirectly further affect the actual flavor dependence of the condensate. We shall come back in more detail 
below on this point in the phenomenological discussion in next section, after establishing our final result for the precise 
condensate values. 


E. Four-loop 0 ( 8 3 ) results 

We finally consider the optimization of the spectral density at four-loop order, the maximal order available at 
present. In fact, the complete standard perturbative expression of our starting expression for the condensate, i.e. the 
next ag order correction to (5.1), is not fully known at present. But it obviously takes the form 

m (qq)QCD P ( m , 9 ) = ^2 w4 ( ) 3 ( C 40im + C 41^m + ^2^ + C 43 im + C44) (5.12) 

where L m = In (in/p) and we choose a convenient overall normalization with respect to lowest order terms in (5.1). 
Now the leading (LL), next-to-leading (NLL) and next-to-next-to-leading (NNLL) logarithms coefficients c 4 o —C 42 are 
easy to derive from RG invariance properties as being fully determined by lowest orders. The NNNLL In m coefficient 
C 43 can also be inferred by RG properties from the available anomalous dimension of the vacuum energy, calculated 
by Chetyrkin and Maier [41], and related to s 3 given in Eq. (5.5). Explicitly, we obtain 

c 40 ^ 4836.74 (4533.33) 
c 4 i ~ -12282.5 (-11292.4) 
c 42 ~ 15606.4 (12648.1) 

c 43 ~ -18588.6 (-15993.5) (5.13) 

where the first and second numbers correspond to n/ = 2 and n/ = 3 respectively. (N.B.: We can obtain the generic 
algebraic values of C 4 *(n/) but these are rather involved and not particularly instructive, so we prefer to keep an 
approximate numerical form for the relevant n/ = 2,3 case in (5.13).) Thus only the nonlogarithmic coefficient C 44 is 
actually unknown at present, and could be quite challenging to compute. But since the nonlogarithmic parts cannot 
contribute to the spectral density, the latter can thus be fully determined at four loops! This gives for the exact 
perturbative four-loop contribution to the spectral density, after taking the logarithmic singularities according to 
Eqs. (2.7),(3.9): 

- 4doop( A ) = ^2 A3 ( c 40 (n f ){2Ll - yT A ) +c 4 i(n/)(^L| - y) + c 42 (n f )L x + y 43 (n/)^ 


(5.14) 
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to be added to the three-loop expression in (5.6). It allows us to calculate the spectral density and the related 
condensate at three successive orders of the variationally modified perturbation, which gives further confidence and 
an important stability and convergence check of our result. 

We obtain at four-loop order once more a unique real common RG and OPT AF compatible solution. (The brute 
optimization results actually give several real solutions for A, as but there are no possible ambiguities since all solutions 
are eliminated from the AF compatibility requirement, except a single one, with as > 0 and L\ < 0 as expected.) 
Explicitly we obtain the optimization results given in the 5tli and 6 th lines of Tables II, III respectively for nf = 2,3, 
where to illustrate the stability the 5th lines correspond to taking the RG Eq. (3.13) at three-loop order, and the 
6 th lines (more consistently) at four-loop order (A being always taken now at four-loop order from Eq. (5.11)). One 
observes a further decrease of the optimal coupling as to more perturbative values, with respect to the three-loop 
results above, as well as the corresponding decrease of L \, meaning that fx is also larger. The stabilization/convergence 
of the results is even clearer for the scale-invariant condensate ( qq)RGi given in the last columns in Tables II, III, 
which at four-loop order has almost no variation upon RG Eq. truncations. 

To better appreciate the very good stability of these results, consider the basic perturbative expression of the 
condensate (5.6) up to four loops in more numerical form (for n/ = 2) and a more standard normalization of the 
coupling: 


~Pqcd( 4 - loop)( A) ^ A 3 (-1 + ^(L A - 0.42) + (^) 2 (9.46 + (29.1 - 25.7 L X )L X ) 
+(^-) 3 (9!.5 + L A (—129 + L x (-288 + 151L A ))) . 


(5.15) 


From this one can easily appreciate that the successive perturbative terms are not small, just like in most perturbative 
QCD series: At successive orders the coefficients grow rapidly (even if partly damped by the decreasing as/ft higher 
powers, provided that as remains rather moderate). In fact for the relevant above values of as — 0.4 — 0.5 roughly, 
and typically L\ ~ —(.7 — . 8 ), depending on the RGOPT order, all successive perturbative terms are roughly of 
the same order of magnitude. Now for the variationally modified perturbation the successive sequences are quite 
different, but before any optimization the resulting series in as has perturbative coefficients that similarly grow at 
successive orders. But the RGOPT mass and coupling optimization manage to stabilize the series, in such a way 
that the discrepancies between the three- and four-loop orders in the final (qq) 1 ^ 3 results are about 2% or less. It 
is important that the optimized sequence has thus clearly further stabilized from three- to four-loop order, to be 
more confident in a quite precise determination of the condensate, although the variation from the lowest nontrivial 
two-loop to three-loop results were already very reasonable, by only ~ 4%. 

It appears that these QCD RGOPT results are more stable than the corresponding ones for the spectral density 
of the yet simpler large -N GN model in Table I. This is a bit surprising a priori , given that direct optimizations, 
not going through the spectral density, give maximal convergence for the large -N GN model [20]. In fact one can 
understand these results as follows. As explained in section IV above, the rather slow convergence for the GN spectral 
density is entirely due to the large and growing factors of 7 r 2p coming from the discontinuities (3.9) at successive 
orders, spoiling the simple form of the series and ‘screening’ the otherwise maximal convergence with the neat so¬ 
lution (4.6). Now although the 7 r 2p coefficients from (3.9) are universal, thus the same for QCD, once combined 
with the original perturbative coefficients of (5.1) their relative contributions with respect to the other perturbative 
terms remain more reasonably of the same order in the QCD case than in the GN case. This is because the original 
perturbative coefficients are comparatively larger in the QCD case. More precisely inspecting the QCD spectral 
density series at three-loop order g 2 in (5.7), the 7 r 2 contribution (last term in the RHS of (5.7)) is roughly twice 
the other non-logarithmic contribution (first term in the RHS of (5.7)). Similarly at the next four-loop g 3 order 
from (5.12) the 7 r 2 contributions are roughly twice the other non-logarithmic contribution (C 43 / 2 ). As long as those 
n 2p contributions remain roughly of the same order of magnitude as the original perturbative coefficients, such that 
some balance can occur from the optimization process, they should produce a moderate disturbance of the observed 
stability. We expect these properties to remain true at even higher orders, because the QCD original perturbation 
coefficients also grow rapidly with the order. 

Lines 4 to 6 of Tables II, III give our direct optimization results for n/ =2,3 and three or four loops respectively. 
To get a final result it could be legitimate to take only the presumably more precise maximal four-loop perturbative 
order available, as is commonly done in most perturbative analysis. However to allow for a more realistic estimate 
of the theoretical error of our results, we will more conservatively consider the difference between the three- and 
four-loop results (but using consistent RG perturbative order in each case) as defining the theoretical uncertainty. 
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VI. EVOLUTION TO HIGHER ENERGY AND PHENOMENOLOGICAL COMPARISON 

In order to get a more precise result it is necessary to take into account the (moderate but not completely negligible) 
running of the condensate values, since the optimal scales obtained are somewhat different at three- and four-loop 
order, though reasonably perturbative, roughly of order fl > 1 GeV. It is anyway necessary to perform a further 
evolution of scale if only to make contact with the more standard scale /i ~ 2 GeV where other (sum rules, lattice, 
etc.) condensate determinations are often conventionally given. 


A. RGOPT (qq)(fj, = 2GeV) results for n/ = 2 and rif = 3 

The procedure to evolve perturbatively the condensate from one to another scale is straightforward since from exact 
RG invariance of m(qq) it is simply given by the inverse of the well-known running of the quark masses, 

/•s(A) ^ ( a \ 

(? 9>(/0 = {qq){n)exp[ / . ( 6 . 1 ) 

JgW P\9) 

Alternatively we may take the values of the scale-invariant condensate (5.10) as obtained in the last columns of 
Tab. II, III and extract from those the condensate at any chosen (perturbative) scale ji' by using again (5.10) now 
taking g = 47 Tas(/x'). This is of course fully equivalent to performing the running from g to // with Eq. (6.1). Since 
all relevant scales jl obtained above are in a fairly perturbative range > 1 GeV, we take a (four-loop) perturbative 
evolution 13 . We choose the highest optimized scales obtained (given by the four-loop results for both n/ = 2 and 
rif = 3 cases) as reference low scale(s): /z re /(2) = 3.45A 2 and /i re /(3) = 3.70A 3 respectively. (N.B.: Given the present 
values of A 3 ~ 340 ± 8 MeV[44], g re f{ 3) happens quite accidentally to be very closely below the charm quark mass 

threshold). For example this running gives a ~ 2% increase of the three-loop {qq)n f= 2 value given in Tab. II and 
quite similarly for rif = 3. Putting all this together we obtain 

-{qq)n f 3 = 2 ^ref( 2) = 3.45A 2 ) = (0.796 - 0.808)A 2 , 

-{qq^Uilirefi 3) = 3.70A 3 ) = (0.773 - 0.796)A 3 , (6.2) 

which are our intermediate results in terms of A and at the respective nf = 2,3 optimal scales, including our esti¬ 
mated theoretical uncertainties, roughly of order 2 %, given by the range of differences between the three-loop results 
(evolved to the scales /x re /) and direct four-loop results. We give results in the form (6.2) in view of possibly more 
precise determinations of A 2j3 in the future. Note that for both n/ = 2,3 the lowest values given in ( 6 . 2 ) correspond 
to the presumably more accurate maximal four-loop results, which gave the g re f values directly from optimisation, 
so without possible extra uncertainties from running. 

Next we perform a final evolution from the low reference scales fi re f{ 2,3) relevant for nf = 2 and rif = 3 respectively 
as given in (6.2), up to the conventional scale // = 2 GeV, where from the present world average A 3 above we find 
bs(2GeV) ~ 0.305 ± 0.004. For n/ = 3 we take into account the four-loop expression of the perturbative matching 
[45] at the crossing of the charm quark threshold. Overall this leads to an increase of the values in (6.2) for |(gq , )| 1 ^ 3 of 
about ~ 4.6% for nf = 2 and 5.3% for rif = 3 (in which taking into accout the charm quark threshold with matching 
relations for as(/r — ?n c ) contributes to ~ —0.3% with respect to a more naive one-step evolution ignoring charm 
threshold effects). More precisely: 

-{qq)]!^ 2 {2GeV) = (0.833 - 0.845)A 2 

-(9<7)n / / 3 = 3 (2GeV) = (0.814 - 0.838)A 3 . (6.3) 

To give a more precise determination for rif = 2 one obstacle is the presently not very precisely known value of the 
basic scale A 2 . In principle it is beyond purely perturbative reach, as it cannot be ‘perturbatively connected’ to the 
more precisely known A 3 value [44]. Our own estimate [22] of A 2 from the pion decay constant gave A 2 — 360^30 MeV, 
while some recent lattice determinations are A 2 — 330 ±45 (staggered Wilson fermions [46]) and A 2 — 331 ±21 (quark 


13 For rif = 2, it is implicitly understood that this evolution is performed in a simplified QCD world where the strange and heavier quarks 
are all infinitely massive, i.e. ‘integrated out’. Otherwise it would not make sense perturbatively to take into account the strange quark 
mass threshold effects on the running. For rif = 3 we can perform a more realistic running, taking into account properly the charm 
quark mass threshold effects on asi^ ~ 77i c ), see below. 
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static potential method [47]). (Incidentally the latter most precise present lattice determination tended to increase 
the central value of A 2 by ~ 15 GeV with respect to former similar determinations [48]). Since lattice uncertainties 
are mostly statistical and systematic, while ours are theoretical errors, it is not obvious how to combine all these in a 
sensible manner. We thus prefer to keep separate estimates. For a representative illustration, combining our present 
results in (6.3) with the above quoted most precise lattice values of A 2 , we obtain 

-<««)^ 3 = 2 (2GeV, lattice A 2 ) ~ 278 ± 2 ± 18 MeV; (6.4) 

where the first quoted error is our intrinsical theoretical error propagated from the one in ( 6 . 2 ), while the second larger 
uncertainty originates from the lattice ones on A 2 . Using instead solely our above quoted RGOPT determination [22] 
of A 2 gives somewhat higher values with larger uncertainties: 

-(&)„%( 2GeV, rgopt A 2 ) ~ 301 ± 2±HMeV. (6.5) 

For rif = 3 the more precisely known A 3 value from many different determinations allows a more precise determination 
of the condensate. Taking the latest world average values [44] as(mz) = 0.1185 ± 0.0006 translating in A|’ a = 340 ±8 
MeV, we obtain 

~{m)nf=z{ 2 GeV,A^ a ) ~ 281 ± 4 ± 7 MeV ( 6 . 6 ) 

where again the first error is our estimated theoretical uncertainty and the second one is from the world averaged A 3 . 
Using alternatively solely our RGOPT determination [22] of A 3 = 317t 2 o MeV, gives instead slightly lower values but 
with larger uncertainties: 


-(QQfJfUi 2 GeV, rgopt A 3 ) ~ 262 ± 4±% MeV. 


(6.7) 


Finally rather than fixing the scale from A, it may be more sensible to give our results for the ratio of the scale- 

1 /3 

invariant condensate with another physical scale, which is a parameter-free prediction. Taking the (gg )results in 
Tab. II, III, and using solely our previous RGOPT results [22] for F/ A 2 and F 0 /A 3 (where F (F 0 ) are the pion decay 
constant for n/ = 2, n/ = 3 respectively in the chiral limit), we obtain: 


-<^> 1/3 


RGI ,rif —2 


F 


= 3.25±0.02t°- 35 


-0.24 


( 6 . 8 ) 


where the first error comes from the present calculation of the condensate while the second ones from taking the most 
conservative range combining linearly three- and four-loop order uncertainty results for F/ A 2 from Eq. (4.28) of [22]. 
A less conservative estimate may be obtained alternatively by taking the range spanned by the maximal available 
four-loop results for F/ A 2 correlated with the four-loop condensate results. This gives 


-\ L i L URGi,n f =2 = 3 26 q. q1 +o 

F 


Similarly for n/ = 3 we obtain: 


-m 


1/3 

RGI ,n f —3 


F 0 


3.04 ±0.04^;^ , 


(6.9) 


( 6 . 10 ) 


where the first theoretical error comes from the condensate with the second one from the conservative range combining 
linearly three- and four-loop order uncertainty on F 0 /A 3 from Eq. (4.30) of [22]. As observed above, the direct results 
from optimization of ~{qq){nf)/h^ lf in Tables II, III show a moderate relative decrease, of about 2 — 3% only on 

{Qg)n f - 3 /A 3 . The effect appears slightly more pronounced, about 7% relative reduction from n/ = 2 to n/ = 3 when 
comparing the central values of ( 6 . 8 ) and (6.10), due to the slight 4% reduction of (central) F 0 relative to F, although 
this result is not clear, being affected with rather large uncertainties. We may finally combine ( 6 . 8 ) and (6.10) to give 


/qqX 1 / 3 t p 

^ f/3 J ’" /=3 - (°- 97 ± °- 01 ) IT ~ (°' 94 ± °- 01 ± °' 12 ) "77 - °- 86 ± °- 01 ± O - 11 ± °- 05 > 

( 99 ) RGI,n f =2 2 


( 6 . 11 ) 
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where all our theoretical errors are combined linearly. In the results on the RHS of (6.11) the first quoted errors are 
the intrinsic RGOPT errors for the present condensate calculation only, and the second larger one is propagated from 
the F/A 2 and F 0 /A 3 RGOPT theoretical errors. We also stress in (6.11) that our results are by construction in the 
strict chiral limit m q 0. The result given for unspecified Fq/F corresponds to the present sole RGOPT condensate 
estimate without extra input from other methods, while the last result uses the present lattice Fq/F estimates [2] 
(with its own uncertainty ~ 0.05 quoted last). 


B. Comparison and discussion 

One may compare (6.4), (6.5) with the latest, presently most precise lattice determination, from the spectral 
density [13] for n/ = 2: (<?<?)„ =2 (a« = 2GeV) = —(261 ± 6 ± 8), where the first error is statistical and the second 
is systematic. Our results in (6.4) are thus compatible within uncertainties, though marginally if taking the other 
RGOPT determination of A 2 in (6.5) tending to be relatively high. Note however that the above quoted lattice value 
from [13] was obtained by fixing the scale with the Fk decay constant rather than using A (moreover with Fk being 
determined in the quenched approximation). It is thus probably more judicious to compare our results for the RG 
invariant ratio (6.8) with theirs [13]: 2.77 ± 0.02 ± 0.04. Overall, recent lattice determinations from various other 
methods lie in the range roughly (qq)l/ 3 4 5 6 7 8 — 2 ~ —(220 — 320) MeV for 71/ = 2 [2], and quite similarly for n/ = 3. The 

most precise lattice n/ = 3 determination we are aware of is (qq)' =3 (2GeV) = — (245 ± 16) MeV [49]. Concerning 
the nf = 3 to n/ = 2 condensate ratio, various lattice results have still rather large uncertainties at present [2] but 
some recent results are compatible with a ratio unity [50]. Our results compare a bit better with the latest ones from 
spectral sum rules [3]: ( uu) 1 2 / 3 ~ — (276 ± 7) MeV (and for the ratio [51, 52]: ( ss)/{uu) = 0. 74 ^ 0 . 12 )- However the 
sum rules results [3] actually determine precisely the current quark masses, extracting then the (uu) value indirectly 
from using the exact GMOR relation (1.1). 

We stress again the rather moderate ?i/ dependence of our result. This is in some tension with the larger estimated 
difference between n/ = 2 and n/ = 3 cases obtained by some authors [18]. Since our results are by construction valid 
in the strict chiral limit, taken at face value they indicate that the possibly larger difference obtained by some other 
determinations is more likely due to the explicit breaking from the large strange quark mass, rather than an intrinsic 
n/ dependence property of the condensate in the exact chiral limit. 


VII. SUMMARY AND CONCLUSION 

We have adapted and applied our RGOPT method on the perturbative expression of the spectral density of the Dirac 
operator, the latter being first obtained from the perturbative logarithmic discontinuities of the quark condensate in 
the MS-scheme. This construction allows successive sequences of optimized nontrivial results in the strict chiral 
limit at two-, three- and four-loop levels. These results exhibit a remarkable stability and empirical convergence. 
The intrinsic theoretical error of the method, taken as the difference between three- and four-loop results, is of order 
2%, while the final condensate value uncertainty is more affected by the present uncertainties on the basic QCD 
scale A, specially with a larger uncertainty for n/ = 2 flavors. The values obtained are rather compatible, within 
uncertainties, with the most recent lattice and sum rules determinations for 71/ = 2 , and our values indicate a moderate 
flavor dependence of the ( qq)n f / A nj ratio. 


[1] M. Gell-Mann, R. J. Oakes and B. Renner, Phys. Rev. 175, 2195 (1968). 

[2] G. Colangelo et al, Eur. Phys. J. C 71 , 1695 (2011); S. Aoki et al., arXiv:1310.8555 [hep-lat]. 

[3] S. Narison, Phys. Lett. B 738, 346 (2014) [arXiv:1401.3689 [hep-ph]]. 

[4] Y. Nambu and G. Jona-Lasinio, Phys. Rev. 122, 345 (1961). 

[5] S. P. Klevansky, Rev. Mod. Phys. 64 , 649 (1992); T. Hatsuda and T. Kunihiro, Phys. Rept. 247 , 221 (1994). 

[6] P. Maris, C. D. Roberts, and P. C. Tandy, Phys. Lett. B 420, 267 (1998); P. Maris and C. D. Roberts, Phys. Rev. C 56, 
3369 (1997). 

[7] K. Langfeld, H. Markum, R. Pullirsch, C. D. Roberts, and S. M. Schmidt, Phys. Rev. C 67, 065206 (2003). 

[8] See e.g. H. G. Dosch and S. Narison, Phys. Lett. B 417 , 173 (1998); M. Jamin, Phys. Lett. B 538, 71 (2002). 



20 


[9] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B147, 385 (1979). 

[10] See e.g. D. Becirevic and V. Lubicz, Phys. Lett. B 600, 83 (2004); P. Hernandez, K. .Jansen, L. Lellouch, and H. Wittig, 
J. High Energy Phys. 0107, 018 (2001); A. Duncan, S. Pcrnice, and J. Yoo, Phys. Rev. D 65, 094509 (2002). 

[11] E. Marinari, G. Parisi, and C. Rebbi, Phys. Rev. Lett. 47, 1795 (1981). 

[12] L. Del Dcbbio, L. Giusti, M. Liischer, R. Pctronzio, and N. Tantalo, J. High Energy Phys. 0602, 011 (2006); L. Giusti, 
and M. Liischer, J. High Energy Phys. 03 (2009) 013; H. Fukaya, S. Aoki, T.W. Chiu, S. Hashimoto, T. Kaneko, J. Noaki, 
T. Onogi, and N. Yamada, Phys. Rev. Lett. 104, 122002 (2010) [Erratum-ifcid. 105, 159901 (2010)]. 

[13] G. P. Engel, L. Giusti, S. Lottini, and R. Sommer, arXiv:1411.6386 [hep-lat]. 

[14] T. Banks and A. Casher, Nucl. Phys. B169, 103 (1980); 

[15] J. Gasser and H. Leutwyler, Annals Phys. 158, 142 (1984); Nucl. Phys. B250, 465 (1985). 

[16] H. Leutwyler and A. V. Smilga, Phys. Rev. D 46 (1992) 5607. 

[17] A. V. Smilga and J. Stern, Phys. Lett. B 318, 531 (1993); K. Zyablyuk, J. High Energy Phys. 0006, 025 (2000). 

[18] S. Descotes-Genon, L. Girlanda, and J. Stern, J. High Energy Phys. 01 2000 041; S. Descotes-Genon and J. Stern, Phys. 
Lett. B 488, 274 (2000); S. Descotes-Genon, N. H. Fuchs, L. Girlanda, and J. Stern, Eur. Phys. J. C 34, 201 (2004); 
V. Bernard, S. Descotes-Genon, and G. Toucas, J. High Energy Phys. 01 (2011) 107; V. Bernard, S. Descotes-Genon, 
and G. Toucas, J. High Energy Phys. 1206, 051 (2012); M. Kolesar and J. Novotny, Nucl. Part. Phys. Proc. 258-259, 90 
[arXiv: 1409.3380 [hep-ph]]. 

[19] J. Bijnens and T. A. Lahde, Phys. Rev. D 71 (2005) 094502. 

[20] J.-L. Kneur and A. Neveu, Phys. Rev. D 81, 125012 (2010). 

[21] J.-L. Kneur and A. Neveu, Phys. Rev. D 85, 014005 (2012). 

[22] J.-L. Kneur and A. Neveu, Phys. Rev. D 88, 074025 (2013). 

[23] C. Arvanitis, F. Geniet, .J. L. Kneur, and A. Neveu, Phys. Lett. B 390, 385 (1997). 

[24] J.-L. Kneur, Phys. Rev. D 57, 2785 (1998). 

[25] V.I. Yukalov, Tlieor. Math. Phys. 28, 652 (1976); W.E. Caswell, Ann. Phys. (N.Y) 123, 153 (1979); I.G. Halliday and P. 
Suranyi, Phys. Lett. B 85, 421 (1979); P. M. Stevenson, Phys. Rev. D 23, 2916 (1981); Nucl. Phys. B203, 472 (1982); J. 
Killinbeck, J. Phys. A 14, 1005 (1981); R.P. Feynman and H. Kleinert, Phys. Rev. A 34, 5080 (1986); A. Okopinska, Phys. 
Rev. D 35, 1835 (1987); A. Duncan and M. Moshe, Phys. Lett. B 215, 352 (1988); H.F. Jones and M. Moshe, Phys. Lett. 
B 234, 492 (1990); A. Neveu, Nucl. Phys. B, Proc. Suppl. B18, 242 (1991); V. Yukalov, J. Math. Phys (N.Y.) 32, 1235 
(1991); S. Gandhi, H.F. Jones and M. Pinto, Nucl. Phys. B359, 429 (1991); C. M. Bender et al., Phys. Rev. D 45, 1248 
(1992); S. Gandhi and M.B. Pinto, Phys. Rev. D 46, 2570 (1992); H. Yamada, Z. Phys. C 59, 67 (1993); K.G. Klimenko, 
Z. Phys. C 60, 677 (1993); A.N. Sissakian, I.L. Solovtsov, and O.P. Solovtsova, Phys. Lett. B 321, 381 (1994); H. Kleinert, 
Phys. Rev. D 57, 2264 (1998); Phys. Lett. B 434, 74 (1998); Phys. Rev. D 60, 085001 (1999); Mod. Phys. Lett. B 17, 
1011 (2003). 

[26] C. Arvanitis, F. Geniet, M. Iacomi, J.-L. Kneur, and A. Neveu, Int. J. Mod. Phys. A 12, 3307 (1997). 

[27] R. Seznec and .J. Zinn-Justin, J. Math. Phys. 20, 1398 (1979); .J.C. Le Guillou and J. Zinn-Justin, Ann. Phys. 147, 57 
(1983); J. Zinn-Justin, arXiv: 1001.0675. 

[28] R. Guida, K. Konishi, and H. Suzuki, Ann. Phys. (N.Y.) 241, 152 (1995); 249, 109 (1996). 

[29] J.-L. Kneur and D. Reynaud, Phys. Rev. D 66, 085020 (2002). 

[30] H. Kleinert, Mod. Phys. Lett. B 17, 1011 (2003); B. Kastening, Phys. Rev. A 68, 061601 (2003); Phys.Rev. A 69, 043613 
(2004). 

[31] J.-L. Kneur, A. Neveu, and M. B. Pinto, Phys. Rev. A 69, 053624 (2004). 

[32] B. Hamprecht and H. Kleinert, Phys. Rev. D 68, 065001 (2003). 

[33] J. L. Kneur, M. B. Pinto, R. O. Ramos, and E. Staudt, Phys. Rev. D 76, 045020 (2007). 

[34] J. L. Kneur, M. B. Pinto and R. O. Ramos, Phys. Rev. C 81, 065205 (2010) 

[35] J. A. M. Vermaseren, S. A. Larin, and T. van Ritbergen, Phys. Lett. B 405, 327 (1997). K. G. Chetyrkin, Nucl. Phys. 
B710, 499 (2005); M. Czakon Nucl. Phys. B710, 485 (2005). 

[36] D. J. Gross and A. Neveu, Phys. Rev. D 10, 3235 (1974). 

[37] Mathematica, version 6, S. Wolfram Company. 

[38] K. G. Chetyrkin and V. P. Spiridonov, Sov. J. Nucl. Phys. 47, 522 (1988); 

[39] J. C. Collins, “Renormalization. An Introduction To Renormalization, The Renormalization Group, And The Operator 
Product Expansion”, (1984) Cambridge University Press, Cambridge, UK. 

[40] K. G. Chetyrkin and J. H. Kuhn, Nucl. Phys. B432, 337 (1994); K. G. Chetyrkin and A. Maier, J. High Energy Phys. 01 
(2010) 092. 

[41] K. G. Chetyrkin and A. Maier, private communication. 

[42] T. Inagaki, D. Kimura and A. Kvinikhidze, Phys. Rev. D 77, 116004 (2008). 

[43] K. G. Chetyrkin, J. H. Kuhn and M. Steinhauser, Comput. Phys. Commun. 133, 43 (2000) 

[44] See the QCD review chapter by S. Bethke, G. Dissertori and G. P. Salam in K. A. Olive et al. [Particle Data Group 
Collaboration], Chin. Phys. C 38, 090001 (2014). 

[45] K.G. Chetyrkin, B.A. Knielil and M. Steinhauser, Phys. Rev. Lett. 79 (1997) 2184 [hep-ph/9706430]; Nucl. Phys. B510, 
61 (1998); K. G. Chetyrkin, J.H Kuhn and C. Sturm, Nucl. Phys. B744, 121 (2006). 

[46] ETM collaboration, B. Blossier et al, Phys. Rev. D 82, 034510 (2010). 

[47] F. Karbstein, A. Peters and M. Wagner, arXiv: 1407.7503 [hep-ph]. 

[48] K. Jansen et al., [ETM Collaboration], J. High Energy Phys. 01 (2012) 025. 

[49] A. Bazavov et al., [arXiv:0910.2966]. 



21 


[50] C. T. H. Davies, C. McNeile, A. Bazavov, R. J. Dowdall, K. Hornbostel, G. P. Lepage and H. Trottier, PoS ConfinementX, 
042 (2012). 

[51] C. A. Dominguez, N. F. Nasrallah, R. Rontsch and K. Schilcher, J. High Energy Phys. 0805, 020 (2008). 

[52] S. Narison and R. Albuquerque, Phys. Lett. B 694, 217 (2010); R. M. Albuquerque, S. Narison, and M. Nielsen, Phys. 
Lett. B 684, 236 (2010). 



